Numerical Study of the Incommensurate Phase in Spin-Peierls Systems 
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We analyze several properties of the lattice solitons in the incommensurate phase of spin-Peierls 
systems using exact diagonalization and quantum Monte Carlo. These systems are modelled by 
an antiferromagnetic Heisenberg chain with nearest and next-nearest neighbor interactions coupled 
to the lattice in the adiabatic approximation. Several relations among features of the solitons and 
magnetic properties of the system have been determined and compared with analytical predictions. 
We have studied in particular the relation between the soliton width and the spin-Peierls gap. 
Although this relation has the form predicted by bosonized field theories, we have found some 
important quantitative differences which could be relevant to describe experimental studies of spin- 
Peierls materials. 

PACS numbers: 74.20.-z, 74.20.Mn, 74.25.Dw 



I. INTRODUCTION 

One-dimensional or quasi-one-dimensional magnetic 
systems show many fascinating properties which con- 
tinue to attract an intense theoretical activity. One of 
these properties is the presence of a spin gap iruantifer- 
romagnetic Heisenberg chains with integer spintl and in 
ladders .□ Another, particularly complex, system which 
presents a spin gap is the spin-Peierls (SP) system. In 
this system a Heisenberg chain coupled to the lattice 
presents an instability at a critical temperature, Tsp, 
below which a dimerized lattice pattern appears and a 
spin gap opens in the excitation spectrum.tj 

The interest in the spin-Peierls phenomena was re- 
cently revived aftei'the first inorganic SP compound, 
CuGe03, was founda This inorganic material allows the 
preparation of better samples than the organic SP com- 
pounds and hence several experimental techniques ca: 
be applied to characterize the properties of this system 
Besides, this compound can be easily doped with mag- 
netic and non-magnetic impurities, leading to a .better 
understanding of its ground state and excitations.Q 

Spin-Peierls systems present also a very rich and inter- 
esting behavior in the presence of an external magnetic 
field. Below the spin-Peierls transition temperature, and 
for magnetic fields H smaller than a critical value H cr (T), 
the system is in its spin-Peierls phase, characterized by 
a gapped nonmagnetic (S z = 0) ground state with a 
dimerized pattern or alternating nearest-neighbor (NN) 
interactions. For T < T tc < T SP , at H = H cr (T) a 
transition occurs from the dimerized phase to a gapless 
incommensurate (IC) state characterized by a finite mag- 
netization, S z > 0. T tc is the temperature of the point 
at which the dimerized, incommensurate and uniform 
phases meet. The dimerized-IC transition was predicted 
by some theoriesQ to be of first order at low temperatures, 
and this is the behavior found in experimental studiesM. 
Other thepries predict that this transition is a second 
order one.E2l 

A simple picture of the dimerized-IC transition can be 



obtained by mapping the Heisenberg spin chain to a spin- 
less fermion system by a Jordan- Wigner transformation. 
The effect of the magnetic field favoring a nonzero S z due 
to the Zeeman energy can be interpreted as a change in 
the band filling of the equivalent spinless fermion system. 
As a result, the momentum of the lattice distortion moves 
away from n as q = (1 — S z /N)ir, where N is the number 
of sites on the chain. However, since umklapp processes 
pin the momentum at i up to a critical field H cr (T), 
the lattice distortion will remain a simple dimerizatiac. 
and the magnetic, ground state will— ['papain a singlet .til 
Theoreticallj'ESEj and experimentali-IEj studies indicate 
that the lattice distortion pattern in the IC phase corre- 
sponds to an array of solitons. A complementary picture 
indicating how a soliton lattice could appear as a conse- 
quence of the finite magnetization in the IC phase is the 
following. Let's assume that the dominant contribution 
to the magnetic ground state comes from a state of NN 
singlets or dimers. An up spin replacing a down spin 
destroys a singlet and gives rise to two domain-walls or 
solitons separating regions of dimerized order which are 
shifted in one lattice spacing with respect to each other. 
Each soliton carries a spin-1/2. Due to the spin-lattice 
coupling it is expected that the lattice solitons are driven 
by these magnetic solitons. 

The soliton formation in spin-Peierls systems has been 
studied analytically by bosonwation techniques applied 
to the spinless fermion model.ES The coupling to the lat- 
tice is treated usually in the adiabatic approximation. 
The resulting field-theory formalism has lead to impor- 
tant results, the most remarkable being the relation be- 
tween—the soliton width and the spin-Peierls gap, £ ~ 
A -1 .£3 Although this formalism has been extended to a 
Heisenberg model with competing NN and next-p£axest- 
neighbor (NNN) antiferromagnetic inter actionai3li3, it 
presents some unsatisfactory features. 

In the fixst place, there are some recent experimen- 
tal resultaij for the soliton width in the IC phase in 
CuGe03 indicating a disagreement with the theoretical 
prediction. Although there might be a contribution to 
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the soliton width coming from magnetic^ or elasticB in- 
terchain couplings which would explain at least partially 
this disagreement, it is also possible that the differences 
could be due to several approximations involved in the 
bosonized field theory. One should take into account that 
these theories are valid in principle in the long wave- 
length limit, and the applicability of their results to real 
materials can not be internally assessed. Then, our first 
motivation to start a numerical study of the 1C phase 
in spin-Peierls systems is to measure the importance of 
these approximations in the analytical approach. 

In the second place, the field theory approach does not 
provide a detailed dependence of the magnitudes involved 
in terms of the original parameters of the microscopical 
models. For example, even for the simplest caseEJ the 
expression obtained for the spin-wave velocity must be 
replaced by the exact one known from Bethe's exact so- 
lution of the Heisenberg chain. In this sense, numerical 
studies could give information about how the relevant 
magnitudes depend on the original parameters without 
further approximations. 

With these motivations, in this article we want to initi- 
ate the study of the incommensurate phase in SP systems 
using numerical methods. These methods give essentially 
exact results for finite clusters, and they can be used to 
check various approximations required by the analytical 
approaches and the validity of their predictions. Besides, 
the numerical simulations provide a detailed information 
of the dominant magnetic and lattice states. In Section 
H we present the model considered and we study several 
features of the soliton formation in the IC phase using the 
Lanczos algorithm. In particular we analyze the effect of 
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NNN interactions on the soliton width. In Section 
we perform Monte Carlo simulations using the world line 
algorithm -which allows us to study larger chains than 
the ones accessible to the Lanczos algorithm- in order to 
reduce finite size effects. 



II. EXACT DIAGONALIZATION STUDY 

The one-dimensional model which contains both the 
antifcrromagnetic Heisenberg interactions and the cou- 
pling to the lattice is: 

N 



H = J y^(l + {u l+ i - m)) Si ■ Sj+i 



i=l 
N 



K 



N 



j 2 ^2 Si ■ s. i+2 + — y^(uj + i - Ui) 



(i) 



where Si are the spin-1/2 operators and u, is the displace- 
ment of magnetic ion i with respect to its equilibrium po- 
sition. Periodic boundary conditions are imposed. The 
first term, which corresponds to the nearest neighbor 
(NN) interactions, contains the spin-lattice coupling in 
the adiabatic approximation. The second term contains 



the AF NNN interactions, which were proposed in Refs. [ 
[j3|^0| to fit the experimental magnetic susceptibility data 
in CuGeC>3 . Several other properties of this maLepsl have 
been reasonably described using this modeLErEl As m 
Ref . [ [n| , we assume for simplicity that the lattice dis- 
tortion docs not affect the second neighbor interactions. 
In principle, the NNN interactions should be corrected 
by a term proportional to (iti+2 — u{) which vanishes in 
the dimerized phase but not necessarily in the incom- 
mensurate phase. This correction should be important 
precisely in the region around a soliton. It is customary 
to introduce the frustration constant a = J2/J. The es- 
timated value of a in CuGeC>3 varies between 0.24 (Ref. 
[ |2(|) and 0.36 (Ref. [ |l9|). In this second case, a is 
larger than the critical value a c « 0.2411 above which in 
the absence, of dimerization a gap opens in the excitation 
spectrum.Ej 

Our purpose is to study numerically Hamiltonian (|l|) 
with exact diagonalization (Lanczos) techniques and by 
Monte Carlo simulations. In this latter case, in order to 
avoid the well-known sign problem due to the frustra- 
tion, we will consider only the diagonal second neighbor 
interaction 



H-2 



N 
i=l 



(2) 



instead of the isotropic NNN interactions (second term 
ofEq. §. 

It is quite apparent that the main numerical difficulty 
is related to the handling of the set of displacements {u{\, 
which in principle can take arbitrary values to describe 
the various distortion patterns present in the dimerized 
and IC phases of the system. These displacements are 
calculated self-consistently by the following iterative pro- 
cedure. First, we introduce the bond distortions defined 
as Si = (ui + \ — Ui). Then, the equilibrium conditions for 
the phononic degrees of freedom: 



d(H) 
d5i 



+ A = 



lead to the set of equations: 



N 



J(S, • S(+i) +KSi- — ■ s m) = 



(3) 



(4) 



which satifies the constraint Y^ - 5i = 0. This constraint 
has been included in Eq. (Q) through the correspond- 
ing Lagrange multiplier A. The expectation values are 
taken with respect to the ground state of the system. 
The iterative procedure starts with an initial distortion 
pattern {i^ ''}, which in general we choose at random. 
At the step n, with a distortion pattern {<$]•" we di- 
agonalize Hamiltonian (|l|) using the Lanczos algorithm 
and compute the correlations (Si ■ Sj+i). We replace 
these correlations in Eq. (0) and the new set {<?■" } is 
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obtained. We repeat this iteration until convergence. Es- 
sentially the same procedure is followed in the quantum 



Monte Carlo algorithm, as it is discussed in Section III 



We have applied this exact diagonalization procedure 
to determine the distortion patterns in the 20 site chain 
at T = 0. In the first place we consider the case of S z = 
0. As mentioned above, this corresponds to a dimerized 
lattice, i.e. Si = (—1) 1 5q. Notice that for this simple 
case, the equilibrium distortion amplitude <5o could be 
determined in an easier way by computing the energies 
of the spin part of Hamiltonian for a set of values of <5o . 
Then, adding the elastic energy and interpolating one 
obtains the minimum total energy We have performed 
this calculation in order to check our iterative algorithm. 

The results for So vs. K, for S* — 0, are shown in 
Fig.l for a = 0.0, 0.2 and 0.4, and Jf = 0.2, and 0.4. It 
can be seen that, as expected, for a > the dimerized 
state is more favorable and this leads to a larger So for a 
given K. To a lesser extent this trend is also present for 

Ji>o. 

The dependence of So with K can be inferred from the 
scaling relation between the energy and the dimerization, 
Eq(So) — Eo(0) ~ Sy^ (plus logarithmic correctionslaa^i. 
v = 2/3, in principle valid for a < a c and small <5o.u&E£l 
Then, it is easy to obtain Sq ~ K~ 3 / 2 , a relation which 
is approximately satisfied by our numerical data. The 
fact that So vanishes at a finite value K of the elastic 
constant, is just a finite size effect. By diagonalizing 
chains of TV = 12, 16 and 20 sites, for a = 0, we have 
verified that K increases with the lattice size, as it can 
be seen in Fig. and it should eventually diverge in the 
bulk limit. 
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1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 
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FIG. 1. Dimerization amplitude vs. elastic constant ob- 
tained by exact diagonalization in the 20 site chain, S z = 0, 
for various values of a = Jij J and Jf • 

Once we have determined the equilibrium distortion 
as a function of K, we are able to compute the singlet- 
triplet spin gap, defined as the following difference of 
ground state energies: 

A = E , dim (S z = 1) - E (S Z = 0) (5) 

It is worth to emphasize that Eo t dim is the ground state 



energy of the system for S z — 1 with the dimerization 
obtained for S z = and the same set of parameters. 
The results of this calculation are shown in Fig. ||. Con- 
sistently with the larger So shown in Fig. 0, the gap in- 
creases with a. The effect of Jf is much weaker than that 
of the isotropic second neighbor interaction which is not 
surprising since the ID ground state magnetic structure, 
with a dominant dimerized state, has essentially a quan- 
tum (off-diagonal) origin. This small increase in A for a 
given K is consistent with the small increase in So shown 
in Fig. |l|. The corresponding scaling relation, A ~ if -1 , 
obtained from the relation between the singlet-triplet gap 



and the dimerization, A 
isfied by our numerical data 

0.4 



c2/3 



is again reasonably sat- 
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FIG. 2. Dimerization amplitude vs. elastic constant ob- 
tained by exact diagonalization for TV = 12, 16, 20 (solid 
squares, diamond and triangles, respectively) and Monte 
Carlo simulations for N — 64 (open dots), with a — 0. The 
inset shows the expected scaling behavior So ~ K~ z ^ 2 for 
TV = 64. 




1/K 

FIG. 3. Singlet-triplet gap vs. elastic constant obtained by 
exact diagonalization in the 20 site chain, for various values 
of a and Jf . The symbols have the same meaning as in Fig. 
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FIG. 4. Dimerization amplitude vs. elastic constant ob- 
tained by exact diagonalization in the 20 site chain, S z = 1, 
for various values of a and J| . 

We now consider the case of S z = 1, which corresponds 
to the incommensurate region just above the dimerized- 
incommensurate transition. We have determined the dis- 
tortion pattern for a 20 site chain using the iterative pro- 
cedure described above. As discussed at the beginning of 
this section, the two solitons or domain walls separating 
dimerized regions are clearly distinguishable. (A typical 
pattern can be seen in Fig. 00 The maximum distortion 
Sq, shown in Fig. ^, presents similar behavior as the one 
shown in Fig. nl corresponding to S z = 0. In particular, 
the fact that So vanishes at a finite K is again due to 
finite size effects. 

In order to compute the soliton width, we use the fol- 
lowing form to fit the numerically obtained distortion 
patterns: 

ft = (-l)^tanh f izipfj tanh ^-*o + ^ ( (6) 

which corresponds to modeling each soliton as an hyper- 
bolic tangent j_as obtained in the analytical approach to 
this problemO The amplitude <5, the soliton width £, 
and the soliton- antisoliton distance d, are the parame- 
ters determined by the numerical fitting. The amplitude 
S should be equal to the maximum distortion defined 
above for well separated solitons, i.e. d ;§> £. The main 
limitation of this calculation arises in the region where, 
for a given a, if is so large that the solitons have a sub- 
stantial overlap in the 20 site chain, and the fitting func- 
tion (^3j) is no longer appropriate. In this case, the elliptic 
sine should be used to describe the soliton lattice. This 
is the region where finite size effects are important, as it 
was discussed above with respect to Figs, [l] and|i|. How- 
ever, this situation is not directly relevant to experiment 
since in real materials the solitons are well-separatedo 
We show in Fig. || the soliton width as a function of 
the gap A for the 20 site chain, for the same values of 



a and J| as before. It can be seen that the there is a 
linear dependence of the soliton width with the inverse of 
the gap. This behavior is consistent with the theoretical 
prediction:E£l 

£ = v s /A, (7) 

where v s is the spin- wave velocity for a < a c . It was 
recently shown that the relation (J?]), originally obtained 
for the unfrusteated chain,E3 is also valid in the presence 
of frustrationJiS For a > a c , A contains a contribution 
from the frustration due to the presence of a gap even in 
the absence of dimerization. 

A linear fitting of these curves in the region £ > 2.5 
gives the slopes 1.87, 1.70 and 1.63 , for a = 0^ 0.2 
and 0.4 respectively. Recently, a numerical studynJ has 
proposed the law: v s — §(1 — 1.12a) in the bulk limit 
for a < a c , From this law one gets v s — 1.57, 1.22, for 
a = 0.0 and 0.2 respectively. We can observe that the 
slopes obtained by fitting the curves shown in Fig. (||) 
are systematically larger than these values of v s . Besides, 
the effect of a is weaker in the numerical data than that 
predicted by Eq. (0). For a — 0.4 > a c w 0.2411, we 
have estimated v s by fitting the excitation dispersion re- 
lation e(k) = E , dim (S z = l,k) - E {S Z = 0,k = 0) 
with the law e(k) 2 = A 2 + v 2 s k 2 + ck 4 around k = and 
5i = 0. For L = 20 we obtained v s — 0.707, a value which 
is also smaller than the slope of the curve £ vs. 1/A for 
a = 0.4 in Fig. (H). This disagreement between the pre- 
diction obtained by the continuum bosonized theory and 
the numerical results could be due to the approximations 
involved in the former or to finite size effects present in 
the latter. The study of much larger lattices than those 
considered in this section will be done in the following 
section using quantum Monte Carlo simulations. On the 
other hand, for the case of J| = 0.4 the slope is actually 
larger (~ 2.1) than the value obtained for the Heisenberg 
chain with NN interactions only. This effect is opposite 
to that of the isotropic NNN interactions and it will be 
further discussed in the next section. 
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FIG. 5. Soliton width vs. singlet-triplet spin gap obtained 
by exact diagonalization in the 20 site chain, S z = 0, for 
various values of a and Jf • 



4 



III. MONTE CARLO SIMULATIONS 

In order to treat longer chains than those consid- 
ered in the Lanczos diagonalization study of the previ- 
ous section, weJiave implemented a world-line Monte 
Carlo algorithms suited to this problem. The parti- 
tion function is re-expressed as a functional integral over 
wordline configurations, where the contribution on each 
imaginary-time slice is given by the product of the two- 
site evolution matrix elements, 

Wi,»+l( T ) = (^i,T^i+l,T | c AtJ,Si S * + 1 | ^ t+ At^Lt+At) 

where Ji = J(l + Si). These matrix elements are the 
Boltzmann weights associated with a bond (i, i + 1) in a 
time step At = 1/ mT in the Trotter direction, where T is 
the temperature and m is the Trotter number. Since the 
exchange couplings depend on the lattice displacements, 
these matrix elements are site dependent. 

We implemented the algorithm with the addition of a 
dynamic minimization of the free energy with respect to 
the lattice displacements. Starting from a given initial 
configuration (random distribution of spins and a dimer- 
ized pattern for the lattice displacements) we typically 
considered 2 x 10 3 sweeps for thermalization. During the 
next 4 x 10 3 sweeps we measured the derivative of the 
magnetic free energy, which, in the limit of T — > 0, is 
given by 

= J«S r S i+1 )) T . (8) 

Leaving 3 sweeps between each measurement for de- 
correlation this produces 10 3 independent values to ob- 
tain the thermal average. With this free-energy gradient 
we corrected the displacements according to (^) and re- 
peated the procedure, including the 2 x 10 3 sweeps for 
thermalization since the spins have to accommodate to 
the new lattice distorsions. Once the displacement pat- 
tern is stabilized within statistical fluctuations — we typ- 
ically considered ~150 iterations, see Fig. ^ — we per- 
formed measurements of several quantities. For this we 
obtained 100 independent groups of 10 3 measurements 
each, following the same procedure as described above, 
i.e., i) thermalization, ii) measurements of d ^ s M and ob- 
servables, and iii) correction of the displacement pattern 
due to statistical fluctuations. 

In our calculations we considered chains of 64 sites 
with periodic boundary conditions and a temperature 
T = 0.05 J. We checked that this value is low enough to 
study ground-state properties by comparison with mea- 
surements at even lower temperatures. On the other 
hand, at higher temperatures the soliton is not observed 
and there is no definite pattern of lattice displacements. 
We took m — 80 for the Trotter number, which is 
large enough to reproduce the Lanczos results on smaller 
chains (see Fig. ||). For some particular quantities like 
the energy gap, which require more precision, we consid- 
ered also to = 160. In addition, comparison with results 



for a longer chain with N = 128 indicates that in the pa- 
rameter range of our calculations the Monte Carlo results 
have no sizeable finite-size effects. 

0.10 i ■ 1 ■ 1 ■ 1 




iteration 

FIG. 6. Minimization of the free energy for a uniform 
dimerized chain. We plot the derivative of the free energy 
and the parameter 8 along the successive iterations. 

In Fig. H we show the Monte Carlo results for the ho- 
mogeneous dimerization of the 64 site chain in the S z = 
subspace as a function of the elastic constant K, together 
with the Lanczos results for smaller chains. Notice that 
in the parameter range considered the 64 site chain does 
not have the finite size effects present for smaller chains, 
namely, the vanishing of Sq for finite values of K. The 
inset shows the expected scaling behavior S oc K~ 3 ^ 2 
discussed in the previous section. As a further check, we 
have also reproduced the scaling behavior of the energy 
gain Eo(Sq) — Eq(0) and gap with <5o with a measured 
exponent v = 2/3 within statistical errors. 
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FIG. 7. Lattice distortion patterns and local magnetiza- 
tions of the 64 site chain obtained by Monte Carlo, for S z = 1 
and different values of K. In every panel the maximum lat- 
tice distortion is normalized to one, so that they cannot be 
directly compared. 
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FIG. 8. Soliton width vs. elastic constant K obtained by 
Monte Carlo simulations in the 64 site chain for Jf = 0.0 
and 0.3, together with Lanczos results for the 20 site chain, 
Jf = 0.0 and 0.4, and a = 0.4. The dashed line corresponds 
to a linear fit to the Monte Carlo results for Jf = 0. 

The soliton structure in the subspace with S z = 1 is 
given in Fig. |t], where we plot the displacement envelope 
Si = (—ySi and the local magnetization (Sf ) , for dif- 
ferent values of the elastic constant K. Notice that the 
displacements are normalized by their maximum values 
(shown in Fig. |J) and the local magnetization by the clas- 
sical value S — 1/2. Consequently, the size of lattice dis- 
tortions in different panels cannot be directly compared. 
For small values of K there is a well defined soliton- 
antisoliton structure in the distortion pattern, with the 
associated local magnetization following a staggered or- 
der. There is a net 1/2 spin density near each domain 
wall, which makes the excess S z = 1. As in the pre- 
vious section, we fitted a two-soliton solution (^|), with 
So = 1 because of the normalization adopted. The re- 
sults for the soliton width £ are shown in Fig. |[ For 
increasing values of K the soliton width grows until the 
displacement profile resembles a sine law (see Fig. [t]). 
This sinusoidal pattern is typical of the soliton lattice, 
observed for large values of S z . It can be seen that the 
scaling £ ~ K obtained in [ |l2| is well reproduced in the 
whole parameter range considered, as indicated by the 
linear fit to the data (dashed line). This figure shows 
that the soliton width for Jf = 0.3 also presents a lin- 
ear dependence with K. These features observed in the 
64 site chain are qualitatively similar to those present in 
the 20 site chain as determined by exact diagonalization. 
Besides, it can seen in this figure that the reduction of £ 
is much stronger when the isotropic NNN is taking into 
account. 

We have performed a simple study on the soliton- 
antisoliton interaction. For this study we fixed the dis- 
tortion pattern to the law (o) with the previously fitted 



value of £, and considered increasing values of d. For small 
K (< 2 J) we found that the total energy becomes a con- 
stant (within statistical fluctuations) when d > 4£, which 
implies that the soliton-antisoliton pairs shown in the left 
panels of Fig. 4 are not interacting. This was confirmed 
by allowing the lattice distortion to evolve starting from 
a pattern like (|J) with an initial separation larger than 
d, which produces the same result for £ and the total 
energy. 
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FIG. 9. Soliton width vs. inverse of singlet-triplet gap ob- 
tained by Monte Carlo simulations for Jf = 0.0 and 0.3. The 
horizontal error bars give the estimated error in the determi- 
nation of the gap A. 

Next, we study the behavior of the soliton width £ with 
the spin-Pcicrls gap A. That is, we compare the quantity 
£ that characterizes the S z — 1 soliton state, with the 
singlet-triplet excitation gap A above the dimerized S z = 
ground state. As shown in Fig. |[ these two quantities 
are inversely related to each other, as discussed in the 
previous section. The slope of the linear fit is 1.9, very 
close to the value 1.87 obtained by exact diagonalization 
of the 20 site chain in the previous section. This result 
confirms the disagreement between the numerical results 
with the analytical prediction pointed out in Section O. 
Also shown in Fig. | are the results for Jf — 0.3. A linear 
fit to these results leads to a slope w 2.3 , i.e. larger than 
the value corresponding to Jf = 0.0. This increase of the 
slope between £ and A" 1 is consistent with the result 
obtained for the 20 site lattice by exact diagonalization 
and Jf = 0.4. This behavior should be contrasted with 
the reduction of the slope found for the isotropic NNN 
interaction. A possible explanation of this behavior could 
be the following. As discussed in the previous section, the 
term Tii^ leads to a smaller increase of the spin gap than 
the fully isotropic NNN interaction. On the other hand, 
the Ising interaction could be more effective in punishing 
the excess {S z ) which appears around a soliton leading 
to a smaller reduction of the soliton width than the one 
caused by the isotropic term, as it can be seen in Fig. ||. 
A more detailed study of the Hamiltonian in the presence 
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of the term of 7i| z is clearly necessary to fully understand 
this behavior. 

Finally, it is possible to estimate the critical value of 
the magnetic field at zero temperature. By adding a 
Zeeman term to the Hamiltonian (Q), —g/j,BS z H (^b- 
Bohr's magneton), H cr may be calculated as: 

H cr = E (S Z = 1)-E Q (S Z =0) (9) 

in units of g^s- Eo(S z = 1) is the ground state energy 
of (|l|), and then H cr < A, which is the value expected of 
a gapped system in the absence of magneto-elastic cou- 
pling. The behavior of H cr as a function of A is shown 
in Fig. [To] for the 64 site chain, a = Jf = 0.0, and for 
the 20 site chain, a — Jf = 0.4. It is apparent a linear 
dependence over all the range studied. .which is in agree- 
ment with the mean-field predictioroO, H cr « 0.84 A. 
However, we obtain a coefficient considerable smaller, 
H cr /A w 0.47, almost independent of a. This value is 
also smaller than twice the soliton formation energy cal- 
culated in Ref. [ [l2|]. The finite value at the origin of 
the curves corresponding to a — Jf = 0.4 is a finite size 
effect. 
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In the first place we have not detected any crossover in 
the behavior of the quantities examined as a, the ratio 
of NNN to NN interactions, becomes greater than a c at 
least for the small chains considered. That is, there are 
only smooth changes as a varies between 0.0 and 0.4. The 
most important effect of the competing NNN interaction 
is a reduction of the soliton width £ as a function of the 
inverse of the singlet-triplet spin gap A. Furthermore, 
the effect of the diagonal term (^) is much less important 
and in some cases even qualitatively different to that of 
the isotropic NNN term. 

Although several functional forms predicted by con- 
tinuum analytical theories have been confirmed by our 
numerical data, there are some important quantitative 
differences. The most important disagreement between 
our numerical results and the analytical predictions is 
related to the coefficient in the relation £ ~ A - , i.e. 
we have obtained a systematically higher value than the 
theoretical value which is the spin-wave velocity. The es- 
timated value of H cr /A is also noticeable smaller than 
the mean-field result and slighly smaller than the pre- 
diction of bosonized field theory. The relevance of these 
numerical results to real SP materials- such as CuGeC>3 
and the recently discovered NaV2C>5,c3 has to be deter- 
mined experimentally. 

The numerical procedures developed in this article 
could be applied to the study of several other proper- 
ties of the incommensurate phase of spin-Peierls systems, 
such as the static magnetization as a function of the 
magnetic field-Xrecently measured in CuGeC>3 by Fagot- 
Revurat et alS3) and the order of the transition from the 
dimerized to the incommensurate phasesJj 
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FIG. 10. H cr vs. spin gap obtained by Monte Carlo in the 
64 site chain for a = J£ = 0.0 and by exact diagonalization 
in the 20 site chain, for a = J% = 0.4. 

IV. CONCLUSIONS 

In this article we have analyzed the magnetic soliton 
lattice in the incommensurate phase of spin-Peierls sys- 
tems using numerical methods. There is a remarkable 
agreement between the results obtained by exact diag- 
onalization using the Lanczos algorithm and those ob- 
tained by quantum Monte Carlo with the world-line al- 
gorithm. The relations among various features of the 
solitons and magnetic properties of the system have 
been determined and compared with analytical results. 
Our starting point is a microscopical model proposed to 
describe several properties of CuGeOa, consisting of a 
ID AF Heisenberg model with nearest and next-nearest 
neighbor interactions. 
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